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Abstract 


A series of analysis methods is proposed to simulate the liquid—gas two-phase and multi-component transport phenomena in the gas diffusion 
layer (GDL) of a proton exchange membrane fuel cell (PEMFC). These methods involve measuring and predicting the two-phase flow properties of 
a GDL, and simulating the two-phase multi-component transport in the GDL. The capillary pressure is measured by the porous diaphragm method 
and predicted by the pore network model. The relative permeability is measured by the steady-state method and predicted by a combination of the 
single-phase and the two-phase lattice Boltzmann method. And the simulation of the liquid—gas two-phase transport is done using the multi-phase 
mixture model. The methods are applied to a carbon-fiber paper GDL to identify the two-phase multi-component transport in the GDL. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Development of proton exchange membrane fuel cells (PEM- 
FCs) has been accelerating for automotive applications since 
the late 1990s. Water management is a crucial factor in improv- 
ing fuel cell performance, which is affected by flooding, gas 
dilution, and membrane dehydration. Reliability of the polymer 
membrane and the catalyst layer is also affected by the presence 
of water [1]. Among water management issues, liquid—gas two- 
phase flow in the gas diffusion layer (GDL) of PEMFCs is an 
important phenomenon because it affects the amount of water in 
the catalyst layer and membrane, and affects reactant transport 
from the gas channels to the catalyst layer. 

Since the GDL of a PEMFC is usually made of electrically 
conductive and opaque material such as carbon fiber, it is hard 
to visualize or measure the liquid condition, phase velocity, and 
species concentration in the GDL. Various simulation efforts 
have been attempted to elucidate the two-phase transport phe- 
nomena in the GDL [2-12] since Wang et al. [2] first introduced 
the concepts of capillary pressure and relative permeability to 
describe two-phase flow in fuel cell GDL assuming Darcy’s law. 
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In general, capillary pressure and relative permeability are 
regarded as the dominant properties of the liquid—gas two-phase 
flow characteristics in a porous medium [13]. Liquid water is 
driven by capillary pressure, which is defined as the difference 
between the gas and liquid phase pressures: 


P; = Pz — P.. (1) 


As the liquid phase pressure changes with the void space occu- 
pied by liquid water, the capillary pressure depends on the liquid 
saturation, s, defined as the liquid volume fraction of the total 
pore space in the porous medium: 


Vi 
S = . 
Vpore 


(2) 


Darcy’s law can be extended for two-phase flow in porous media 
as described in the following equation: 


AP, (3) 


where k, is the relative permeability of each phase defined as the 
ratio of the permeability of the phase at a given saturation level 
to the absolute permeability of the porous medium. 

Although a fundamental understanding of two-phase flow ina 
GDL is of great importance in analyzing the effect of the flood- 
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Nomenclature 


molar concentration (mol m~?) 
capillary number 

mass diffusivity of species (m? s~!) 
lattice velocity vector 

distribution function 

total interaction force 

interaction force between fluid—fluid 
interaction force between fluid—solid 
body force 

body force per unit mass 

interactive strength between phases k and k 
interactive strength between k phase and solid 
phase 

Green’s function 

mass flux of liquid phase (kg m~? s7!) 
current density (A m~?) 

phase 

relative permeability 

permeability (m7) 

length (m) 

molecular weight (kg mol~!) 

mass flux (kg m~? s7!) 

index for solid wall 

connectivity number 

pressure (Pa) 

pore curvature (m) 

gas constant, 8.134 J mol`! K7! 
Reynolds number 

saturation of liquid water 

time 

temperature (K, °C) 

velocity (m s7!) or non-dimensional 
common velocity 

volume (m?) 

location vector 

neighboring site of x 

in direction x, y, Z 


Greek symbols 


B 


LQErY OR 


a Qab 


net water molar flux per proton molar flux through 
a PEM 

convection corrector factor 

porosity 

mobility of a phase 

viscosity (kg m7! s7!) 

contact angle (°) 

kinematic viscosity (m? s7!) 

density (kg m~?) 

surface tension force (N m7!) or non-dimensional 
relaxation time 


Subscripts 


a 
c 


unknown value 
capillary 


eff effective 

eq equilibrium 

g gas phase 

i direction 

Int intrinsic 

l liquid phase 
pore pore in a porous medium 
s solid phase 

sat saturated value 
0 reference state 
Superscripts 

H20 water 

N2 nitrogen 

O2 oxygen 

a species 


ing phenomenon on the PEMFC performance, only a limited 
number of studies on two-phase flow property measurements 
for a GDL used in a PEMFC have been published [14,15]. As 
an alternative approach, the properties have been empirically 
derived from lithological experiments [1,13,16-18]. The simu- 
lations presented in Refs. [2—12] have been attempted with those 
properties from porous media other than an actual GDL as stated 
above. It is therefore desirable to identify a method that can elu- 
cidate the phenomena in GDL resulting from the effects of such 
properties as the fibrous structure and wettability. 


2. Measurement and prediction methods 


This present study uses a carbon-fiber paper as the GDL 
material that has been treated with a 5 wt.% solution of poly- 
tetra-fluoro-ethylene (PTFE) because a GDL is commonly a 
porous medium whose porosity is high enough for gas trans- 
port, with highly hydrophobic wettability, which allows for easy 
drainage of liquid water. Fundamental properties of the GDL 
were obtained by general measurement methods: the contact 
angle, porosity, pore size distribution, and absolute permeabil- 
ity were, respectively, measured by the Wilhelmy or droplet 
method, pycnometer method, mercury porosimetry, and capil- 
lary flow porosimetry [19]. 

The following subsections describe a series of analysis meth- 
ods to elucidate the transport phenomena in a GDL. These 
methods involve measuring the two-phase flow properties of a 
GDL, predicting those properties, and simulating the two-phase 
multi-species transport in a GDL. 


2.1. Capillary pressure measurement 


Fig. 1 shows a schematic diagram of a device for measur- 
ing the capillary pressure versus saturation (P.—s) curve that 
is based on the porous diaphragm method [13,20]. A Toray 
TGP-H type carbon-fiber GDL disk sample is placed inside 
the chamber, and sealing materials are applied on top and bot- 
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Fig. 1. Schematic diagram of P.—s curve measurement device. 


tom surfaces of the sample near the perimeter in order it does 
not interfere with liquid imbibition. Liquid water is forced into 
the porous sample from below by air pressure at room tem- 
perature. The measurement is started by monitoring an abrupt 
change of the liquid water pressure when we observe, thorough 
the transparent sample holder, air bubbles in the liquid water 
below the sample disappear. The capillary pressure and satura- 
tion are, respectively, measured from the liquid water pressure 
and amount of water contained in the pore space of the sample. 
The measurement is done after it is confirmed that the value of 
the liquid pressure has stabilized within a range that the value 
does not affect the P.—s curve. The P.—s curve is determined by 
repeating this liquid water imbibition procedure in small incre- 
ments at each state until liquid water penetrates the GDL. The 
gas phase pressure in the calculation of capillary pressure is 
assumed to be the atmospheric pressure. The saturation is deter- 
mined by precisely measuring the mass of liquid water forced 
into the GDL sample and calculating the ratio of the liquid water 
volume converted from this liquid water mass to the total pore 
volume. 


2.2. Capillary pressure prediction 


We applied pore network model [21,22] to predict the P.—s 
curve of the GDL for two reasons. Firstly, the approach does 
not require high computational cost because the model idealizes 
the pore morphology and topology as a pore network consist- 
ing of pores and throats. Secondly, the GDL properties, such 
as wettability, pore size distribution, and pore connectivity, are 
easily modified as parameters in the calculation. This approach 
is applicable because the process of a phase intrusion can be 
regarded as quasi-static at each saturation value. On the other 
hand, the relative permeability calculation, described in the latter 
section, requires consider the flow in an actual complex struc- 
ture to estimate the pressure drop. The pore network can be 
constructed based on the pore size distribution of the GDL as 
measured by mercury porosimetry (Fig. 2) and the connectivity 
number. The connectivity number is extracted by applying an 
image processing technique that employs a thinning algorithm 
[23,24] to the image that has been obtained by the microfocal 
X-ray CT of the carbon-fiber paper GDL. The resolution of the 
image is | wm. Fig. 3 shows the image of the Toray TGP-H-060 
GDL, which consists of non-woven carbon fibers each with a 
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Fig. 2. Pore size distribution of GDL (Toray TGP-H-060). 


diameter of about 7 um, and the respective porosity is 0.79. The 
x-axis in the figure shows the thickness direction of the GDL, 
which we believe to be the dominant direction in terms of trans- 
port between the gas channel and the catalyst layer. As seen 
in the figure, the fibers are stacked in layers in the x direction, 
and no particular structural difference can be seen in the y and 
z directions. Thus, we assumed that the GDL shows different 
permeability in the x direction, but that it has the same per- 
meabilities in the y and z directions. For the above reason, we 
have measured or predicted absolute and relative permeability 
in the x direction. Fig. 4 shows the result of applying the thin- 
ning algorithm. The gray objects show the fibrous structure of 
the GDL, and the distributed dots show the skeleton of the pore 
network, the flow paths in the GDL. The connectivity number is 
calculated by counting the average number of skeleton lines at 
each connection point. The pore network consists of 120 pores 
in each of 3 directions (120 x 120 x 120). And the diameter of 
each pore is determined by random number based on the pore 
size distribution. The throats are also generated by random num- 


« 


Fig. 3. Microfocal X-ray CT image of carbon paper (Toray TGP-H-060). 
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Fig. 4. Skeleton of pore network in microfocal X-ray CT image of GDL (Toray 
TGP-H-060). Analyzed region 120 pm? in size. 


ber in order that the averaged number of throats corresponds to 
the connectivity number. 

Applying the Young—Laplace law as in Eq. (4) to the pore 
network enables the P.—s curve to be calculated: 


20 COS be 
P= ——_., (4) 


r 


The contact angle, pore curvature based on the pore size dis- 
tribution, and surface tension of water are used in this equation. 
Water is forced from the inlet to the outlet of this pore network 
according to Eq. (4). Calculation parameters are described in 
Table 1. The contact angle is the advancing contact angle mea- 
sured by the Wilhelmy method. We chose this method because 
we assumed that the advancing contact angle is appropriate for 
this liquid intrusion process. 

As an example, Fig. 5 shows the water intrusion behavior in 
the pore network with contact angle 0. = 162° at s=0.091 and 
P. =6805 Pa. It can be seen that there are some points where 
liquid water has been deeply infiltrates, which may lead to liquid 
water breakthrough at low saturation level. 


2.3. Relative permeability measurement 


Various measurements of relative permeability have been 
conducted by the steady-state test method for reservoir rock, 
sand stone, and other porous media [25-27]. We also used the 
steady-state method to measure the relative permeability of the 
gas phase. A schematic drawing of the test apparatus is shown 
in Fig. 6. A Toray TGP-H type carbon-fiber GDL test piece is 


Table 1 

Calculation parameters for pore network model 
Contact angle, 6. (°) 162 

Surface tension, o (N m7!) 70.25 x 1073 
Connectivity number, N 4.9 


Fig. 5. Water distribution in GDL in applied pore network model (6, = 162°). 


sandwiched between similar GDLs on the inlet and outlet sides. 
The GDL on the inlet side ensures homogeneous distribution of 
liquid water in the planar direction, while the one on the out- 
let side minimizes the effect of the outflow boundary. Liquid 
is injected first and then a constant flow rate of air is estab- 
lished at room temperature. The pressure difference is measured 
after the pressure is stabilized. The GDL test piece is removed 
immediately after the pressure measurement so that the weight 
can be measured to deduce the saturation. The measured pres- 
sure difference, air flow rate, and saturation are combined as 
in Eq. (3) with the absolute permeability measured by capil- 
lary flow porosimetry to derive the relative permeability of the 
gas phase. The relative permeability of the gas phase versus 
saturation (krg—s) curve can be obtained in this manner. 


2.4. Relative permeability prediction 


The two-phase lattice Boltzmann method (TLBM) can be a 
useful tool for studying the complex behavior of two-phase flow, 


Fig. 6. Schematic diagram of relative permeability measurement device. 
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Table 2 
Calculation conditions for SLBM used in Steps 2, 4, and 5 


Table 4 
Calculation conditions for TLBM used in Step 3 


Property Value 


Property Value (in lattice units) 


Velocity model 
Calculation domain, Lx x Ly x Lz (um) 


D3Q15 model 
122 x 122 x 122 


Density of gas, pg (kg m7?) 1.0 
Density of liquid, pı (kg m~?) 1000 
Viscosity of gas, 4g (Pas) 21 x 1076 


370 x 1076 
7.0 x 1074 


Viscosity of liquid, u (Pas) 
Inlet velocity, uy (ms~!) 


such as coalescence and the breakup of water in a porous medium 
with the Young—Laplace law [28]. Although it is possible to sim- 
ulate the arbitrary contact angle and capillary number of a water 
droplet by applying TLBM to GDL, it is difficult to represent 
the real properties of liquid water and air at the same time in 
complex geometry [29]. Such properties include the viscosity, 
density, and surface tension. To resolve this issue, we propose 
a method that combines TLBM described in Refs. [28-31] and 
the single-phase lattice Boltzmann method (SLBM) described in 
Ref. [32] to predict the relative permeabilities of both phases ver- 
sus the saturation of GDL (A;g—s, k-s). The method is described 
by the following steps: 


e Step 1: reconstructing a voxel image of a GDL from micro- 
focal X-ray slice images. 

e Step 2: calculating the absolute permeability of the voxel 
image by the SLBM approach. Details of the conditions are 
shown in Table 2. Inlet velocity is specified for this step. 

e Step 3: calculating the two-phase flow by the TLBM 
approach. The governing equations for the TLBM are tab- 
ulated in Table 3 and the Green’s functions are adjusted to 
have specified contact angles beforehand. Details of the con- 
ditions are shown in Table 4. With these parameters, the order 
of capillary number Ca based on the velocity and viscosity of 
the liquid phase is less than 1075 and Reynolds numbers for 
both phases Re), and Reg based on the mean pore diameter 
are less than 1077. These are set in order that surface tension 
force is more predominant than viscous force and the flow can 
be treated as Darcy’s flow. These assumptions are based on 
rough estimation of those non-dimensional numbers at actual 
fuel cell operating condition. A body force is applied as the 


Table 3 
Governing equations for TLBM 


Velocity model 
Calculation domain, Lx x Ly x Lz 


D3Q15 model 
122 x 122 x 122 


Density of gas, pg 80 

Density of liquid, p; 160 

Viscosity of gas, Hg 16.7 (relaxation time of gas, 
Tg =1.15) 

Viscosity of liquid, 4 50.1 (relaxation time of liquid, 
T= 1.44) 


Contact angle, 6. 135° (interaction strength 
regarding wettability, 

8gs = —8ls = —0.025) 

18.92 (interaction strength 
regarding the surface tension, 
8g = 0.001) 

Body force, gy 0.001 


Surface tension, o 


driving force and the periodic boundary condition is applied. 
The water distribution in the GDL can then be obtained at a 
certain saturation level. 

e Step 4: in calculating the permeability for the water distri- 
bution obtained in Step 3 by the SLBM approach, liquid 
water and fibers are regarded as solid in this step. The rel- 
ative permeability of the gas phase (krg) can then be obtained 
by calculating the ratio between this permeability and the 
absolute permeability calculated in Step 1. 

e Step 5: performing similar calculations for the relative per- 
meability of the liquid phase (kn). The gas phase and fibers 
are regarded as solid in this step. 


Fig. 7 shows the fibrous structure and water distribution in 
GDL when liquid saturation s=0.216 in Step 3. For the reason 
stated in Section 2.2, we have measured or predicted absolute 
and relative permeability in the x direction. We assume that the 
computational domain size is examined. It is because the poros- 
ity of the domain matches to the macroscopic averaged value of 
the whole domain, and the absolute permeability calculated by 
SLBM on this computational domain also agreed with the value 
estimated from manufacturer’s data. 


Equations Value (k=1, g) 


FE) (OM KD 


Hatert D- fin = 


Lattice Boltzmann equation 


Tk 


Macroscopic velocity pu = pkU' + | Fx, where u’ 


Fluid/fluid interaction 


Fluid/solid interaction 


x 


Body force F3: = Pkg 


Fu = -0Y Gesla, xn — x), where ns = { 


PRUK/ Tk K(x, te; 
2 u = J fee , and Fy = Fig + Fox + F 3, 
i 


SS pelt , 


Fu = p10) Y Ga, xoa — x), where Gig (x, x’) = 
x! k 


Ekk Ix—x|=1 
gy/V3 |x —x'| = V3 
0 otherwise 
i Sks Ix -—x|=1 
1 solid 
0 fluid and Gxs(x, x’) = 8ks / V3 Ix —x’| = /3 
0 otherwise 
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Fig. 7. Fibrous structure and water distribution in GDL using TLBM. s=0.216, 
120 pm? computational region. It flows along with x direction. 


2.5. Transport simulation 


Our two-phase flow simulation model is based on the multi- 
phase mixture model developed by Wang and Cheng [16]. It 
solves a set of macroscopic equations for multi-phase multi- 
component flow in a porous medium and can simulate the 
two-phase flow and species transport in a GDL. In this study, 
the equations are arranged for a steady-state isothermal immis- 
cible two-phase system involving a liquid water phase and gas 
phase, which consists of oxygen, water vapor, and nitrogen as 
shown in Table 5. Physical values without subscript, which are 
density, velocity, concentration, pressure, and viscosity, are for 
the mixture of the gas and liquid phases as defined in the table. 
The velocities in the formulation are superficial velocities in a 
porous medium. Gravity is neglected since pore size is so small, 
as shown in Fig. 2, that the capillary effect is dominant. And 
the effective gas phase diffusion coefficient is modified by the 
porosity as given in the following equation: 


pp = eo De (5) 


Table 5 


Flux equivalent of Saturated air 


current and transport with liquid 
from membrane 
n f 
Ce 
g 
n H9 Cathode GDL 
C H,O 
O, 8g 
x=0 x=L 


Interface with 
gas channel 


Interface with 
catalyst layer 


Fig. 8. 1D computational region and boundary conditions of transport simula- 
tion. Mass fluxes applied to continuity equation, water mixture equation, and 
oxygen equation at x=0. Gas phase pressure, oxygen concentration, gas phase 
water concentration, and saturation are specified at x = L. Computational region 
divided into 100 control volumes. 


The set of equations in Table 5 employs Darcy’s law for each 
phase as described in Eq. (3), and assumes interfacial thermal 
and chemical equilibrium between the phases. The saturation, 
which is calculated by the liquid water concentration, affects the 
capillary pressure, and the relative permeability. The equations 
can be arranged in a simple mathematical form that is suitable 
for numerical simulation and are implemented in a 3D numer- 
ical solver of advection diffusion equations based on the finite 
volume method described in Ref. [33]. 

A steady-state 1D computation was carried out for a cath- 
ode GDL with boundary conditions for water generation and 
oxygen consumption corresponding to current density J at one 
end (x=0) and fixed pressure, fixed species concentration, and 
fixed liquid saturation at the other end (x = L) as shown in Fig. 8. 
The physical properties, dimensions, and boundary condition 
values are summarized in Table 6. The capillary pressure and rel- 
ative permeability characteristics (Pe—s, krg—s, and k,—s) of the 
Toray TGP-H type GDL were applied as empirical correlation 
with saturation and will be described in Section 3. The absolute 
permeability is measured by the capillary flow porosimetry. Val- 
ues for the oxygen diffusion coefficient, liquid water viscosity, 
liquid water density, and gas viscosity were taken from the liter- 
ature [4]. The computational region is divided into 100 equally 


Governing equations for liquid—gas two-phase mixture with oxygen and nitrogen in gas phase and phase changing water 


Governing equation (phase k=1, g) 


— CHO 


PR" pı 
Continuity V - (pu) = 0, where p = > PkSk, Pg = 3 Cz M*, pu = > PkUk, S1 = mo Some and cm0 = iho 
sat 
=1 
1 f1 u krk 2 krk 
Momentum — | -V . (puu) | = V - (pet Vu) — VP — —u, where u = pv, v = — ,VP= AKV Pk, and Ap = 
e le K Vk i 
cr Cg 
Species a=02, H20: V - (yauC®*) = V - [eD A s)YCg]- V. 7 D. jı |. where C* = ` Ch sk, 
1 g 
p 1 Cy Pàg H20 Pl ; ÀgA1 Pegi 
20 = 4, | A p , = ——=—, cih0 = , Co#HO — 0, j = K ~— —* Vs, and P, = Py — P 
YH20 Ca G F Ya#H20 Pl- s) mmo’ “1 Ji v as 1 c g 1 
a=N2: Cy? = Pe = 


RT j 
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Table 6 

Physical properties, dimensions, and boundary conditions 

Property Value 
GDL thickness, L (um) 200 
GDL porosity, £ 0.8 
GDL absolute permeability, K (m7) 9x 107! 
Saturated water vapor concentration, C#2° (mol m7) 16.1 
Temperature, T (°C) 80 

Net water transport coefficient, 8 0.5 

Gas pressure at x= L, Pg|x=1 (Pa) 202,650 
Oxygen concentration at x=L, ce lr=z (mol m~?) 11.2 
Saturation at x= L, s|,=1 0.2 
Current density, J (A m~?) 20,000 


spaced control volumes. The net water flux transferred through 
the membrane from anode to cathode is assumed to be constant 
and the water mass flux at x=0 yields 


mo_ (1 IM™9 
n? = (3+4) E (6) 


It should be noted that variables that are solved here are the 
mixture of the gas and liquid phases and boundary conditions are 
applied as the appropriate mixture value. In case of the pressure 
boundary condition, the liquid phase pressure at the gas channel 
end (x= L) is higher than the gas phase pressure by the capillary 
pressure as in Eq. (7): 


Pilx=_ = Peles = Pe(s|x=L). (7) 


This can be implemented by applying an appropriate mixture 
pressure P, which can be evaluated by the definition in Table 5, 
as the boundary condition. Note also that the combined mass 
flux of water and oxygen is applied to the continuity equation at 
x=0. 


3. Results and discussion 
3.1. Capillary pressure 


Fig. 9 shows a comparison of the measured results for the 
Toray TGP-H type carbon-fiber paper GDL with the calculated 
results. The experiments were repeatedly performed on samples 
that had been treated with 5 wt.% solution of PTFE. The repeated 
experiments showed similar P,—s curves, though some variation 
occurred due to the difficulty in measuring the small amount of 
liquid saturation. A typical result is shown in the figure. The data 
is plotted only to s =0.1. This is because, during the experiments, 
we observed liquid droplets penetrated to the top surface of the 
sample within the liquid saturation range from 0.1 and 0.2. And it 
could not be clearly determined by observation at what saturation 
the breakthrough occurred. The calculation was performed using 
a pore network model with the GDL properties. 

It can be seen from the comparison that the calculated data 
are similar to the measured data in magnitude, shape and break- 
through point of the P.—s curve. The curve near the breakthrough 
point has a gentle slope, indicating that the liquid water could 
easily intrude when the capillary pressure would reach near the 
breakthrough point. Since the majority of pores in the GDL are 


Capillary Pressure, P, [Pa] 


0) 0.05 0.1 0.15 
Liquid Saturation, s 


Fig. 9. Comparison of measured capillary pressure with calculated data. 


concentrated in the vicinity of the most frequent pore diameter 
from the mercury porosimetry data (Fig. 2), this corresponds to 
the breakthrough capillary pressure. 

The calculation method was therefore successful in express- 
ing a P.—s curve using the GDL properties, thus leading to a 
physical interpretation of the relationships between the GDL 
properties and the curve. Furthermore, the simulation method is 
able to predict the P.—s curve for GDLs of various properties, 
such as wettability, pore size distribution, and pore connectivity. 


3.2. Relative permeability 


A comparison between the calculated relative permeability 
(krg—s and k,—-s) and experimental relative permeability data 
(krg—s) is presented in Fig. 10. A Toray TGP-H type GDL treated 
with a5 wt.% solution of PTFE was used for this experiment and 
calculation. In contrast with capillary pressure prediction, static 
contact angle 0. = 135°, measured by the droplet method, was 
used in this simulation (Table 4) because dynamic contact angles 


x) g 
B O Kr_gas (Calc.) 
o 08 - | O kr_liquid(Cale.) 
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A a EA ®, RO gap 
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Fig. 10. Comparison of measured relative permeability with calculated data. 
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Fig. 11. Cross-sectional images of water distribution in GDL: (a) hydrophobic 
case (6, = 135°) and (b) hydrophilic case (0. =45°). 


such as advancing and receding contact angles are spontaneously 
represented as a numerical artifact by the TLBM approach. The 
accuracy of these contact angles is not examined at this point. 
The relative permeability of the liquid phase near s =0, and that 
of the gas phase near s = 1, could not be predicted by this simu- 
lation; this was because the paths of the liquid phase in the GDL 
were not connected at a low level of saturation from the inlet to 
outlet, and the paths of the gas phase were not connected at a 
high level of saturation from the inlet to outlet. In other words, 
there is a possibility that non-connectivity of the liquid phase 
observed in TLBM leads to negligible relative permeability in 
the low saturation range. 

The predicted gas phase relative permeability correctly cap- 
tured the shape of the curve of the measured results and showed 
a similar magnitude to those at a low level of saturation. It 
is obvious that the absolute permeability of GDL can be esti- 
mated with SLBM correctly. And it can be assumed that the 
path of the gas phase in GDL is formed correctly because the 
relative permeability of the gas phase is agreed with experi- 
ment. In other words, the liquid path in GDL is well represented 
as well. We therefore assumed that the prediction method can 
be applied to the relative permeability of the liquid phase 
successfully. 

To see the effect of wettability on the liquid formation in 
detail, TLBM computation was performed for different con- 
tact angles. Fig. 11 shows the effect of wettability on the liquid 
water distribution in the GDL; Fig. 1 1(a) shows a cross-sectional 
image of the water distribution in a hydrophobic GDL with con- 
tact angle 6e = 135° (already shown in Fig. 7), while Fig. 11(b) 
shows this in a hypothetical hydrophilic GDL with contact angle 
0: =45° at the same saturation level. It can be seen that liq- 
uid water tended to be present in large pores as spheres in 
the hydrophobic GDL, and in small pores as thin films in the 
hydrophilic GDL. This can be explained as follows. When the 
gas phase pressure is assumed to be constant, the liquid phase 
pressure in a small pore is lower than that of a large pore from 
Eqs. (1) and (4). The low liquid phase pressure in the small 
pore draws liquid from big pores. On the contrary, the liquid 
phase pressure is higher in a small pore for the hydrophobic 
case. And it pushes out the liquid phase to big pores. In other 
words, the wetting phase, which is the liquid phase in hydrophilic 
case and the gas phase in hydrophobic case, tends to stay in 
small pores. 
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Fig. 12. Measured capillary pressure and correlated model. 


3.3. Transport simulation 


The empirical equations for capillary pressure and relative 
permeability were formulated by correlating the measured P,—s 
data from Fig. 9 with the measured k,g—s data from Fig. 10, and 
with the predicted k,;—s data also from Fig. 10. These correlated 
models were then applied to a numerical computation. It should 
be noted that the liquid saturation is defined in the measurement 
and the predictions as the ratio of the total liquid water volume 
to the total pore volume in a GDL sample although liquid water 
may be unevenly distributed along the normal direction to the 
liquid injection surface of the sample. On the other hand, the 
liquid saturation in the transport simulation is considered to be 
the ratio of the average local volume of liquid water to the aver- 
age local volume of pore. It is also important to note that the 
capillary pressure could not be measured or predicted above a 
liquid saturation level of about 0.1, and the relative permeability 
of the liquid phase could not be predicted under a liquid satura- 
tion level of 0.2, as described in the previous section. In order for 
these parameters to be applied to the macroscopic approach, we 
have extrapolated the data in Figs. 9 and 10 to fit the empirical 
equations. The fitted correlations are shown in Figs. 12 and 13 
and they are described as Eqs. (8)—(10). The correlated capillary 
pressure is as follows: 


0<s <0.05 


0.05<s<sy<1 


8 
—10*s — 6x 10° (8) 


E x 106s? — 2.5 x 10°s 
P: = 

where s, is the upper bound of the applicable range of Eq. (8). 
The capillary pressure curve decreases linearly after s> 0.05, 
there has to be inflection point somewhere before s= 1 and the 
capillary pressure has to decrease sharply. This is because it 
needs big liquid pressure to force liquid into the remained small 
pores near s=1. The inflection point could not be measured 
experimentally in this study and the value s, is not clear at this 
point. And the correlated relative permeabilties are as follows: 


kı = MAX(1.089(s!> — 1) + 1, 0), (9) 


kg =(—s). (10) 
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Fig. 13. Measured and predicted relative permeability and correlated model. 


We define MAX(A, B) to denote the greater of A and B in Eq. 
(9). And the relative permeability of the liquid phase is 0 when 
saturation is small in this expression. As shown in Table 6, the 
liquid saturation at the gas channel end (x= L) was set to be 0.2. 
Since we found no direct in situ measurement on saturation in the 
GDL at the channel interface in past research, a realistic value 
was set for the first try to demonstrate the simulation method. 
This value was set because an in situ observation of a PEMFC 
shows that the liquid droplets are covering some portion of the 
GDL surface [34]. 

Fig. 14 shows computational results for the species concen- 
tration and liquid saturation. The water vapor concentration is 
constant because isothermal condition is assumed. The liquid 
saturation is high at the catalyst layer end (x=0) due to water 
generation, while the oxygen concentration is low due to con- 
sumption. Similar states can be seen in the pressure profiles of 
the both liquid and gas phases as shown in Fig. 15. The liquid 
phase pressure is high, due to the water generation, while the 
gas phase pressure shows an opposite profile, due to the oxygen 
consumption. The pressure in the figure is the pressure differ- 
ence from the reference state, which is the gas phase pressure at 
the gas channel end (Pg|,=1). Note that the liquid phase pressure 
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Fig. 14. Distributions of species concentration and liquid saturation. 
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Fig. 15. Pressure distributions of liquid phase and gas phase. 


is higher than the gas phase pressure by the capillary pressure 
as in Eq. (7). 

The mass flux and intrinsic local mass average velocity [35] 
were extracted from the computational result at the middle point 
(x=L/2) to identify in detail the transport phenomenon in the 
GDL. The intrinsic velocity of phase k and species a in the gas 
phase are deduced by Eqs. (11) and (12), respectively: 


Uk 
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ans ay 12 
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The results are summarized in Table 7, which shows that the 
gas phase mass flux is negative and that it flows from the gas 
channel end to the catalyst layer end, as opposed to the liquid 
phase, which flows from the catalyst layer end to the gas chan- 
nel end. Oxygen, due to its larger concentration gradient, moves 
much faster than the gas phase. Water vapor moves slowly to the 
catalyst layer along with the advection of the gas phase mixture 
which is caused by the gas phase pressure gradient although 
no diffusion due to concentration gradient occurs. Although 
nitrogen does not move in the macroscopic sense, the diffu- 
sion velocity, which is the difference between the velocity of 
the species and the gas phase velocity, is positive. This means 
that the positive diffusion and negative advection are balanced. 
Compared to the liquid phase water mass flux, transported gas 
phase flux is around half in the reverse direction. And the gas 
phase water flux is negligible. 


Table 7 
Mass flux and velocity in GDL at x= L/2 


Phase Mass flux (kg m~? s~!) Intrinsic velocity (ms~!) 
Gas phase 
Mixture Pgüg = —1.67 x 10-3 UInt,g = — 1.44 x 10-3 
O2 n92 = —1.66 x 107° UQ? a = —7.73 x 107° 
H20 nt = —1.57 x 10-5 umo = —8.44 x 1075 
N2 nw? =0 ny, = 0 
Liquid phase 
Water piu =3.75 x 10-3 Unt =2.40 x 1075 


136 T. Koido et al. / Journal of Power Sources 175 (2008) 127-136 


Although this transport simulation model of the two-phase 
flow properties of GDL needs to be developed further to incor- 
porate the thermal effect and the effect when coupled with other 
layers, it does provide insight into the hidden phenomena of 
a GDL. Experimental validation of the calculated saturation is 
also desirable. 


4. Conclusions 


We have proposed a series of methods to analyze the transport 
phenomena involved in liquid—gas two-phase multi-component 
flow in the GDL of a PEMFC. These methods were applied to 
measure the capillary pressure and relative permeability of the 
gas phase in an actual carbon-fiber paper GDL, to predict the 
capillary pressure and relative permeability for both phases in 
the GDL, and to simulate two-phase multi-component flow in 
the GDL. The methods are applied to a carbon-fiber paper GDL 
considering the microstructure and the wettability of the GDL. 

The predicted capillary pressure and relative permeability 
correctly estimated the shape of the curve of the measured 
results and showed similar magnitude. Transport simulation was 
carried out with the properties to identify the two-phase multi- 
component transport in the GDL. 
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